Explicit radius for nearly parallel vortex filaments in equilibrium 1 
Geophysical research has focused on flows, such as ocean currents, as two dimensional. 

Two dimensional point or blob vortex models have the advantage of having a Hamilto- 
nian, whereas 3D vortex filament or tube systems do not necessarily have one, although 
they do have action functionals. On the other hand, certain classes of 3D vortex mod- 
els called nearly parallel vortex filament models do have a Hamiltonian and are more 
accurate descriptions of geophysical and atmospheric flows than purely 2D models, es- 
pecially at smaller scales. In these "quasi-2D" models we replace 2D point vortices with 
vortex filaments that are very straight and nearly parallel but have Brownian variations 
along their lengths due to local self-induction. When very straight, quasi-2D filaments 
are expected to have virtually the same planar density distributions as 2D models. An 
open problem is when quasi-2D model statistics behave differently than those of the re- 
lated 2D system and how this diff erence is manifested . In this paper we study the nearly 
parallel vortex filament model of 



Klein et al 



( 199a ) in statistical equilibrium. We are 



able to obtain a free-energy functional for the system in a non-extensive thermodynamic 
limit that is a function of the mean square vortex position R 2 and solve explicitly for R 2 . 
Such an explicit formula has never been obtained for a non-2D model. We compare the 
results of our formula to a 2-D formula of iLim and Assadl (|2005h and show qualitatively 



different behavior even when we di sallow vortex br aiding. We further confirm our results 



Klein et al 



using Path Integral Monte Carlo (|Ceperlevl (|l995)) without permutations and that the 



(1995|) model's asymptotic assumptions are valid for parameters where these 



deviations occur. 
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1. Introduction 

In the study of geophysical turbulence, statistical equilibrium theories have focused 
mainly on 2-D models, which are known to have significant differences from the realistic 
3-D. Of particular interest in the 2-D regime is the vortex density distribution in the 
plane and how this density carries over to a quasi-2D regime. In the strictly 2D setting, 
Lim and Assad have variationally derived a low-temperature formula for the variance in 
vortex position, R 2 , for the Onsager model in the Gibbs canonical ensemble with equal 
strength vortices in which 

where N is the number of filaments all with unity circulation, (3 is the inverse temperature, 
and ji is the Lagrange multiplier for the angular momentum. They have confirmed this 
formula and shown that the density is essentially uniform and axisymmetric using Monte 
Carlo simulations indic ating that R 2 is the only significant moment at low-temperature 



( Lim and Assac 



(|2005l )). While we assume that extremely straight filaments have a den- 
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Figure 1. Shown here in top-down projection, nearly parallel vortex filaments, at low density 
and high strength of interaction, are well-order ed into a 2 - D tri angular lattice known as the 



Abrikosov lattice from type-II superconductors IjAbrikosovl |l957|)). This figure shows how the 
quasi-2D model is essentially a 2-D model for these parameters. 

sity profile similar to that of the 2-D model and in turn a similar R 2 (Figure [1} , we 
must also accept that at some magnitude of curviness the R 2 for the nearly parallel 
model differs significantly from that of the fully parallel model. The entropy internal of 
each filament due to self-induction alone must at some point affect the profile. What the 
difference in profiles is and for what parameters it is significant are open problems. 

In this paper we address the problem of when R 2 of the 2-D model agrees with that 
of the nearly parallel vortex filament model of Klein, Majda, and Damodaran (KMD) 
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Klein et al.l (|l995l )). the simpliest model of 3-D interacting filaments. We solve explicitly 

for R 2 in the case of nearly parallel vortex filaments using a simple mean-field approxima- 
tion for the interaction potential and a spherical constraint. We show that with current 
mathematical knowledge a spherical constraint is crucial to solving for R 2 , making ours a 
novel application of the constraint and the first analytical expression for R 2 in a non-2D 
model of vorticity. 



In a previous paper ([Andersen and Liml (|2006l )) we solved for R 2 in a broken segment 



model. Here we address the continuum model and find that our mean-field energy func- 
tional is equivalent to that of a quantum harmonic oscillator with a spherical constraint 
and a constant force. Therefore, to solve for R 2 , we a pply a modified version of the sp her- 
ical quantum methods of Hartman and Weichman (jHartman and Weichmanl (|1995I )) to 



derive the free-energy of the system / as a function of R 2 and minimize / w.r.t. R 2 , 



giving 



R z 



!3 2 aN + y/f3 4 a 2 N 2 + 32a/^ 



(1.2) 



where a is the core-strength parameter and affects how curvy the lines are (jKlein et al 



(1995)), N is the number of filaments, and fj, and are the same as in Equation 11.11 
We show that equation 11.21 in the a — > oo limit (making the lines perfectly straight), 
approaches equation ll.il ffl3")l. 

To confirm our quasi-2D formula and that we do not violate the KMD model's asymp- 
totic assumptions (vortices very straight and f ar apart enough so they do not braid), we 



perform Path Integral Monte Carlo (PIMC) (jCeperlevI (| 19951) ) on the KMD system in 
the Gibbs canonical ensemble at different inverse temperatures, j3. We adapt PIMC to 
our m odel by eliminating the permut ation sampling used fo r quan tum bosons (jCeperlev 



( 19951 )) and type-II superconductors (jNordborg and Blatterl (| 19981 )) and adapt it to sam 



pie exactly the normally distributed, non-interacting filament case of self-induction plus 
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angular momentum so that the Monte Carlo's rejections depend only on interaction 

energy 

For our chosen parameters, the Monte Carlo shows that the model is valid down to 
a certain value of [3 beyond which vortices begin to braid, and it shows that filament 
straightness assumptions hold for even smaller (3. Further, it confirms that (3 values do 
exist where no braiding occurs and our quasi-2D formula for R 2 (Equation II. 2| gives a 
better prediction than Equation ll.il a novel result (Q. 



2. Problem and Model 



Geophysical vortex models are almost exclusively two dimensional because, on a plan- 
etary scale or even the scale of several hundred or thousand kilometers, oceanic and 
atmospheric currents have only a tiny vertical component. Models like the 2-D Onsager 
N-Point- Vortex Gas model are nice and simple to treat mathematically. 

The Onsager Gas Model uses the Hamiltonian, 



h n° = -^A i A J -log|V i -^| 2 , 



(2.1) 



where and Xj are the circulation constants for point vortices i and j and ipi and ipj 
are their positions in the complex plane. The Gibbs canonical ensemble, 

exp(-f3H N - fiI N ) 



Gn — 



where 



and 



Z N = dip 



Z N 



dipN exp (-/3H N - fj,I N ) 



N 



(2.2) 



(2.3) 



i 

represents the statistical distribution of the vortices in the gas with an inverse tempera- 



ture f3 and angular momentum chemical potential /x (| Onsager 



19491 )). From this model, 
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Figure 2. Shown here in perspective projection, nearly parallel vortex filaments are quasi-2D. 
Lim and Assad were able to derive variationally a low-temperature (large [3) formula for 
the second moment of the density R 2 (Equation II. ip and describe the shape of the den- 
sity profile as be i ng ne arly uniform with a sharp density cutoff at 2R when Xi = AjVz, j 



(Lim and Assac 



(2005)). Because the Onsager model is a Hamiltonian system such an 
approach is possible. 

On the other hand, three dimensional models usually have no Hamiltonian repre- 
sentation and models need to simulate the Biot-Savart interaction directly and obtain 
stationary states from the Euler equations, a computationally intensive and analytically 
difficult (if not intractable) task. Despite this difficulty, 3-D models are important to the 
study of realistic geophysics. 

At a certain scale vortex filaments become not quite 3-D but not quite 2-D either and 
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what we term a quasi-2D model is applicable, one that takes into account the variation 

in the vorticity field in the vertical direction but does not allow for hairpins or loops or 

tangles in filaments (Figure [2]) ■ 



A restricted quasi-2D model due to Klein, Majda, and Damodaran (jKlein et al 



dl995h ) 



derived from the Navier-Stokes equations, represents vorticity as a bundle of N filaments 



that are nearly parallel to the z-axis. They have a Hamiltonian, 

„l N - r. , , s 2 pL n n 



J ° k=l J ° k=l i>k 



(2.5) 



where i/jj(o~) = x j(°~) + is the position of vortex j at position a along its length, 

the circulation constant is same for all vortices and set to 1, and a is the core structure 



constant (jKlein et al 



(|1995I )V The position in the complex plane, ipj(o~), is assumed to be 
periodic in a with period L. Self-stretching is neglected in the model although stretching 
due to interaction occurs. Because this system does have a Hamiltonian, we can form 
a Gibbs distribution (Equation 12. 2|) for it and study it in statistical equilibrium with 
conservation of angular momentum, 

In = J2 do-\M°-)\ 2 - (2-6) 
i J o 

This allows us to study the density profile of the vortices. 

For a large range of parameters, the probability distribution for this quasi-2D model 
shows very little difference from that of the fully 2-D model of Onsager because the 
variation in the z-axis direction is too small to have any effect. Therefore, the model begs 
the question of when the density profile begins to behave differently, when do 3-D effects 
come into play, when does this model succeed while the 2-D model fails. 

We hypothesize that there exists and set out the find parameters (a,f3,/j,) such that the 
most important statistical moment, the second one i? 2 , differs in the two systems. (The 
mean is always the origin in this system so by default the second is the most significant 
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moment.) Furthermore, we hypothesize that for at least part of this range of parameters 
the assumptions of the quasi-2D model hold, where these assumptions state that 

(a) for a given rise along a filament e, the variation in the complex plane is order e 2 , 
where e<l 

(6) vortices do not form braids 

(c) the core-size 5 <C e, which we assume. 

We begin by deriving a mean-field, spherically constrained, explicit expression for R^ D , 
Equation 11.21 and we show that, in the a — > oo limit, causing the vortices to become 
perfectly straight and parallel, R^ D — > R\ D: that is that as the quasi-2D model becomes 
fully 2D the expression for R^ D becomes 2D as well (32]). We then verify our expression 
and the second hypothesis, that the a s sump tions are not broken, using Quantum Monte 



Carlo methods of Ceperley (|Ceperlevl l| 19951 )) ($3} 



3. Mean-field Theory 

The partition function for the quasi-2D system, 

Z N = f Di> x ■ ■ ■ [ D^ N exp (S N ) , (3.1) 



where Dtpk represents functional integration over all Feynman paths and Sn = —(3Hn — 
/j,In is the action functional gives us the free energy for the most probable macrostate: 

F = --\ogZ N . (3.2) 

P 

This free energy for the quasi-2D Hamiltonian (Equation 12. 5p cannot be found analyt- 
ically with current mathematical knowledge. In fact, the partition function for the 2-D 
point-vortex system cannot be found for a finite number of vortices. The reason for the 
difficulty comes from the interaction term (second term in Equation l2.5|) . While the other 
terms, the self-induction (first term in Eq. 12. 5p and the conservation of angular momen- 
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turn term, are quadratic and yield a normally distributed Gibbs distribution that we can 

functionally integrate, the logarithmic term must be approximated. 

Although Lions and Majda have made such an approximation and rigorously derived 

a mean-field PDE for th e probability distribution of the vortices in the complex plane 



( Lions and Maidal (|2000h ). their PDE takes the form of a non-linear Schroedinger equa- 
tion that is not analytically solvable, again because of the interaction term. While a 
discussion of their derivation is beyond the scope of this paper, their existence proof 
justifies our applying a much simpler mean-field approximation to interaction. 

For very straight filaments that arc well-ordered, we can take the mean-field action to 



be 



S 



mf 
N 



da 



' N 

E 

.k=i 



da 



dip k {o-) 



Da 



NB 



N 



N 



£l0g|V> fe (<T)| 2 -/i^|^(cT)| 



k=l 



k=l 



(3.3) 



where the second term in the integrand is now mean-field and the other two are the same 
as before. This mean-field app r oxima tion derives from the work of Lim and Assad for 



point- vortices (jLim and Assadl (|2005j )). Intuitively it says that on average the distance 
between vortices is also the distance of a vortex from the center. Based on Monte Carlo 
simulations, this approximation works extremely well even for high levels of filament 
variation. 

Before we begin to calculate the partition function, we must deal with another problem: 
we still cannot integrate Equation 13.11 using this action because the interaction term 
(second term in I3.3() is still a function of ijj and still inside the integral, so we make 
another approximation, adding a microcanonical constraint on the angular momentum, 



5 



da [\4,{a) 



R 2 



(3.4) 
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that has integral representation 
8 



-it [\ip(a)\ 2 - R 2 



(3.5) 



This approximation only applies when fi is large and fluctuations in the angular momen- 
tum (and consequently R 2 ) are fairly small, but it allows us to remove the interaction 
term from the functional integral 13. II 

Now the spherical-mean-field partition function looks like 



Zp ' = J D ^i---J D^ N exp(S N ) ^-exp£ da-ir[\i,(a)\ 2 -R 2 ] 



(3.6) 



jj and, since the exponents are all positive definite, we can interchange the integrals and 
combine exponents, 



7 smf 
J N 



where the action functional is 



dr 
2^ 



5' 



smf 
N 



D-ipjy exp 



/ qsmf 



(3.7) 



N 



fe=i 



(3.8) 



and 



Sk = 



PLN\og{R 2 )/A- 



daa[3 



Da 



+ {it + 2 A1 )|^( ( j)| 2 - iR 2 T {a) 



(3.9) 



is the single filament action. Because tpk are statistically independent for all k, we let 



S = S'fcVfc and S^ 1 = NS, which makes the partition function, 



7 smf 
J N 



N 



(3.10) 



f The Hartman and Weichman paper (H artman and Weichmari 



1 1995)) confines the wave- 



function to the surface of a sphere using an infinite number of Dirac delta functions and a 
functional integral over r, which here would be like confining ip(cr) to a cylinder of radius R. 
We do not apply their method in this way but confine the mean of tp, L _1 f Q daip(a), to that 
cylinder with only one Dirac delta. 
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We need to have the partition function in the form Z s ^ a ^ = j dre~ A T f to use steepest 

descent (Appendix [Bj , and so we define the non-dimensional free energy as a function of 

ir, f[ir] = f3F, 



f\ir] 



log 



Dtp exp (S) 



which allows us to re-write the partition function as 



ysmf _ 
N ~ 



dr_ 

2^ 



exp(-7V/[ir]) 



(3.11) 



(3.12) 



This form of the partition function we can solve with steepest descent methods, but first 
we need a formula for /. 

We can evaluate /[«t], the energy of a 2-D quantum harmonic oscillator with a c onsta nt 



force, with Green's function methods but do not go into it in this paper (jBrownl (|1992 )). 
Let us make a change of variables A = ir + 2/Lt and the non-extensive scaling (3' = (3N and 
a' = a/N, which will become necessary when we take the limit TV — > oo. After evaluating 
the integral, Equation l3.11i (see Appendix |C|) . the free-energy reads 



/[A] =Ln- \h\R 2 - (3'L log(i? 2 )/4 - In 6 " 

2 ( e -wi _ \y 



(3.13) 



where u> = \J\/(a! [3') is the harmonic oscillator frequency. 

Now that we have a formula for / we can apply the saddle point or steepest descent 
method. (For discu ssion of this me t hod s ee Appendix [5] as well as the original paper 
of Berlin and Kac ([Berlin and Kad (|1952l )').) The intuition is that, as N — > oo in the 
partition function, only the minimum energy will contribute to the integral, i.e. at infinite 
N, the exponential behaves like a Dirac delta function, so 



foo = lim InZV^f = fin], 



(3.14) 



where n is such that df[X]/dX\ n = (jHartman and Weichmanl (|1995l ) 



Berlin and Kac 



dl952j) ). 
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First we can make a simplification by ridding Equation 13.131 of R 2 . We know that R 2 

will minimize / and so df jdR 2 — 0. Since 

>>f = {ji-\/2)L-L/3'/(4R 2 ), (3.15) 



dR 2 



R2 = W^WY (3 ' 16) 



Substituting the left side of 13. 161 for R 2 in Equation l3.131 we get 



/[A]=/3'L/4+,/-A 7 i + 21og 



— l0g 4(^A/2) - (3 ' 17) 



a'/3' 

We could take the derivative of Equation I3.17I and set it equal to zero to obtain r\. 
However, doing so yields a transcendental equation that needs to be solved numerically. 
Since our goal is to obtain an explicit formula, we choose to study the system as L —> oo. 
In fact such an approach is justified by the assumptions of the model that L have larger 
order than the rest of the system's dimensions. (If this were a quantum system, this 
procedure would be equilvalent to finding the energy of the ground state. Hence, we call 
this energy f grn d ) Taking the limit on Equation I3.17I yields the free energy per unit 
length in which rj can be solved for 

f 9 rnM = ^ + Vv/WP 7 ) ~ ^ log L./„/o^ . ( 3 - 18 ) 



where 



of which we take 



4 v " v ' ' 4 ° \A(p-r]/2) 



r) = 2(i- l(3'(~f]' 2 a' ± V/?' 4 a' 2 + 32a'/3», (3.19) 
8 



rj = 2/i - \(3'{-(3' 2 a' + V/3' 4 a' 2 + 32a'/3» (3.20) 
8 

as giving physical results. 
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With 77 explicit, we can give a full formula for R 2 , 



R 2 = 



-f3' 2 a' + y/jJ'^a' 2 + 32a'/?' /i 

_ I3' 2 a' + V/3' 4 a' 2 + 32a'/?> 
8a'/3'/x ■ 

Through several approximations, we have obtained an explicit formula for the free 
energy of the system and R 2 . Before we move on to the numerical verification, we can 
ask a few things about it. One question is how 13.211 compares to the formula that Li m 
and Assad derived variationally for point vortices, Equation ll.il (|Lim and Assadl (|2005)). 
We repeat that formula here: 

R\ D = (3.22) 

Because the quasi-2D system has only the self-induction term to differentiate its behav- 
ior from that of the Onsager Vortex Gas, if we make the lines perfectly straight through 
some limit, then the statistics of the quasi-2D system should be identical to that of the 
2-D gas. In the quasi-2D system a' controls the straightness of the lines. The larger a' 
the more inclined the lines are to be straight. Therefore, if we take a' — > 00 limit on 
Equation 13.211 we should get the 2-D formula [3.221 

R 2 (a' 00) = lim 4 



-j3'2 a i + y/pia 12 + 32a' /3'n 

= hm , 

7^0 -pri + ^/3'4 + 32 7 ^' M 

where 7 = a' -1 . We use L'Hopital's rule taking derivatives of the top and bottom 



R 2 (a' -» 00) = lim 

k ' 7-0 if/9'*-, .,_ . ,., . ,_ 

V NR 

(3.23) 



±(/3' 4 + 327/3'm)"3 32/3'/x 
A/3' 2 (3' N/3 



16/3' n 4^ 4^ 

which agrees with Equation l3.22l That these equations agree indicates that, like the point 
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vortex equation. 13.211 is a low temperature formula but can predict behavior for all levels 
of internal filament fluctuation. 

There are qualitative differences in the two equations as well since they are different 
for finite a'. The 2-D equation. 13.221 is linear in (3, but Equation 13.211 is not. In fact, 
a simple calculation (not shown here) tells us that, while 13.221 decreases linearly with 
decreasing (3', Equation 13. 211 begins to increase at a (3' = f3' where 

(3' 3 = (3.24) 
a 

Therefore, for changing (3, the 2-D Equation 13.221 is a straight line, and the quasi-2D 
Equation 13 . 2 1 1 has a "v"-shape, decreasing with f3' to (3' and then increasing again. Since 
the free energy in Equation l3.18l is smooth for all (3 > provided a > 0, [i > 0, this "v"- 
shape does not indicate a phase transition but rather that at some (3 the internal filament 
fluctuations overtake the conflict between the angular momentum and interaction terms 
as the major contributor to the value of R 2 . 

We can also measure quantitative differences between 13.221 and 13.211 We calculate the 
error as a function of the parameters: 

E = - * 2D , (3.25) 
which, after some simple algebra, gives the (3 value at which a particular error appears, 

A a'E(E + l) y ' 

This equation is the most useful for numerical verification because it tells us exactly 
which parameters give how much error E, and error is what we aim to detect since that 
will tell us when 3-D effects become noticeable. In the next section we address whether 
we can detect the error in the Monte Carlo and without violating the assumptions of the 
model. 
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4. Monte Carlo Comparison 

We apply Monte Carlo in this paper to the original quasi-2D model with Hamiltonian 
12.51 to verify two hypotheses: 

(a) that the 3-D effects, namely the Equation I3.21[ predicted in the mean- field are 
correct 

(b) that these effects can be considered physical in the sense that the model's asymp- 
totic assumptions of straightness and non-braiding are not violated. 

We are caught between three requirements: to make the vortices wavy enough to show 
the effects we predict but not so wavy as to violate either the assumptions of the model 
or the assumptions we made in our derivation in $3] The assumptions of the model 
must be violated at some small (3 because decreasing [3 increases fluctuation and vortex 
density, meaning that vortices will begin to braid and ultimately become too wavy to be 
considered asymptotically straight. (There are no other constraints on the model besides 
the filaments having very small core-size, which we assume.) Furthermore, because our 
assumptions make Equation 13.211 a low-temperature prediction, we need entropy-save 
that of internal fluctuations-to be small enough to neglect. In this section, we show not 
only that such a regime exists but that the mean- field equation, Equation I3.21[ predicts 
the mean square vortex position in the Monte Carlo, which we call B? MC , to better 
accuracy than the 2-D equation. 

Research on flux-lines in type-II superconductors has yielded a close correspondence 
between the behavior of vortex filame nts in 3-space and paths of quantum bosons in 



(2+l)-D (2-space in imaginary time) (jNordborg and Blatterl (|1998f ) 



Senetal 



( 20011) '!. 



This work is not related to ours fundamentally because type-II superconductor flux-lines 
do not have the same boundary conditions. They use periodic boundaries in all directions 
with an interaction cut-off distance while we use no boundary conditions and no cut-off. 
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Besides the boundaries, they also allow flux-lines to permute like bosons, switching the 
top endpoints, which we do not allow for our vortices because it would create unacceptable 
tangling. However, despite the boundary differences, the London free-energy functional 
for interacting flux- lines is closely related to our Hamiltonian 12.51 and so we can apply 
Path Integral Monte Carlo (PIMC) in the same way as it has been applied to flux-lines. 
(For a discussion of PIMC and how we apply it see Appendix |A"1) 

We simulated a collection of N = 20 vortices each with a piecewise linear representation 
with M — 1024 segments and ran the system to equilibration, determined by the settling 
of the mean and variance of the total energy. We ran the system for 20 logarithmically 
spaced values of (5 between 0.001 and 1 plus two points, 10 and 100. We set a — 10 7 , 
[i = 2000, and L = 10. We calculate several arithmetic averages: the mean square vortex 
position, 

N M 

Rl IC = {MN)-^Y.\^ k )\ 2 i 

2 = 1 fc=l 

where k is the segment index corresponding to discrete values of a, the mean square 
amplitude, 

N M 

A 2 = (ilOO-^EW*)- Vi(l)| 2 , (4.2) 

i=l fe=l 

which measures how wide the vortices are in top-down projection, the mean square 
amplitude per segment, 

N M 

a 2 = (MN)- 1 E \Mk) ~ Mk + 1)| 2 , (4-3) 
i=l fe=i 

where ipi(M + 1) = tpi(l), and the mean square nearest neighbor distance, 

JV 




Measures of Equation 14 . 1 1 correspond well to Equation 1 3 ■ 2 1 1 in Figure [3] whereas Equa- 
tion [2221 continues to decline when the others curve with decreasing f3 values, suggesting 
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that the 3-D effects are not only real in the Monte Carlo but that the mean-field is a 

good approximation with these parameters. 

In order to be considered straight enough, we need 

L 10 



a < 



(4.5) 



M 1024 

Non-braiding requires that d 2 > A 2 . According to Figure braiding occurs around 
(3 < 0.038 and straightness holds for all (3 values, shown in Figure We do not need to 
show that these conditions hold for every instance in the Monte Carlo, only close to the 
average, because small probability events have little effect on the statistics. We would 
also like to point out that although we cannot guarantee that the statistical results for 
(3 £ (0.002, 0.038) are valid, we cannot dismiss them as being unphysical. Whether the 
statistics for (3 < 0.038 are physical is a subject for future research. 

In Figure S]we show a close-up of the most relevant [3 values in the range (0.038, 0.16) 
where there is greater correspondence between the Monte Carlo results and the quasi-2D 
formula in Equation 13 . 2 1 1 than with Equation ll.il These results are most likely to make 
a good physical prediction. 



5. Related Work 



As mentioned in the previous section, simulations of flux lines in type-II su perconduc- 



tors u s ing the PIMC me thod have been done, generating the Abrikosov lattice (jNordborg and Blatter 



<ll998h 



Sen et al 



(2001)). However, the superconductor model has periodic boundary con- 



ditions in the xy-plane, is a different problem altogether, and is n ot applicable to trapped 



Klein et al 



( 1995) have been done to date 



Kevlahanl |2005J) 



fluids. No Monte Carlo studies of the model of 

and dynamical simulations have been confined to a handful of vortices 
added a white noise term to the KMD Hamiltonian, Equation 12. 5[ to study vortex re- 
connection in comparison to direct Navier-Stokes, but he confined his simulations to two 
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Figure 3. The mean square vortex position, defined in Equation [4TTJ compared with Equations 
I3.21I and II. II shows how 3-D effects come into play around /3 = 0.16. That the 2D formula 
continues to decrease while the Monte Carlo and the quasi-2D formula curve upwards with 
decreasing (5 suggests that the internal variations of the vortex lines have a significant effect on 
the probability distribution of vortices. 

vortices. Direct Navier-Stokes simulations of a large number of vortices are beyond our 
c omputational capacitie s . 



Tsubota et al 



( 20031 ) has done some excellent simulations of vortex tangles in He- 
4 with rotation, boundary walls, and ad hoc vortex reconnections to study disorder in 
rotating superfluid turbulence. Because vortex tangles are extremely curved, they applied 
the full Biot-Savart law to calculate the motion of the filaments in time. Their study did 
not include any sort of comparison to 2-D models because for most of the simulation 



vortices were far too tangled. The inclusion of rigid boundary walls, although correct for 
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Figure 4. The mean square vortex position, defined in Equation [4TTJ compared with Equations 
I3.21l and ll.ll in close-up view between f3 — 0.16 and j3 = 0.038 where we can guarantee validity. 

the study of He-4, also makes the results only tangentially applicable to the KMD system 
we use. 

Our use of the spherical model is recent and has also been applied to the statistical 
mechanics of macroscopic fluid flows in order to obtain exact solutions for quasi-2D 



turbulence (|Liml (|2006), 



Lim and Nebusl |200J) ) 



Other related work on the statistical mechanics of turbulence in 3-D vortex lines 
can be found in 



Lions and Mai da 



Fland oli and Gubinellil (|2002f ) and 



Berdichevskvl (|1998f ) in addition to 
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Figure 5. This figure shows that the mean amplitude per segment (Equation l4.3[l . meaning the 
distance in the complex plane between adjacent points on the same piecewise-linear filament, 
is much less than the segment height L/M — 10/1024, indicating that straightness constraints 
hold for all j3 values. 

6. Conclusion 

We have developed an explicit mean-field formula for the most significant statistical 
moment for the quasi-2D model of nearly parallel vortex filaments and shown that in 
Monte Carlo simulations this formula agrees well while the related 2-D formula fails at 
higher temperatures. We have also shown that our predictions do not violate the model's 
asymptotic assumptions for a range of inverse temperatures. Therefore, we conclude that 
these results are likely physical. Wc have avoided braided vortices. However, the model 
will admit braiding as long as vortices are sufficiently straight and far apart enough. 
In future we will address the problem of vortex reconnections which can change the 
qualitative behavior of the system for some parameters. 
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Figure 6. The mean square amplitude (Equation I4.2[l becomes greater than the mean square 
nearest-neighbor distance (Equation I4.4|l at about /3 = 0.026 indicating that vortices begin to 
braid at this point. Because braiding statistics are not covered in our quasi- 2D model, we cannot 
guarantee that the results for f3 < 0.038 are physical. At this point braiding statistics are the 
subject of future research. 

Appendix A. Path Integral Monte Carlo method 

Path Integral Monte Carlo methods emerged from th e pa t h inte gral formulation in- 



vented by Dirac that Richard Feynman later expanded (|Zed (|2003l )). in which particles 
are conceived to follow all paths through space. One of Feynman's great contributions 
to the quantum many-body probl em was the mapping of path i ntegrals onto a classical 



system of interacting "polymers" (jFeynman and Wheeled (|1948l )). D. M. Ceperley used 



Feynman's convenient piecewise linear formulation to develop his PIMC method which 
he successfully applied to He-4, generatin g the well - know n lambda transition for the first 



time in a microscopic particle simulation (jCeperlevi (jl995t) V Because it describes a system 
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of interacting polymers, PIMC applies to classical systems that have a "polymer" -type 
description like nearly parallel vortex filaments. 

PIMC has several advantages. It is a continuum Monte Carlo algorithm, relying on 
no spatial lattice. Only time (length in the z-direction in the case of vortex filaments) is 
discretized, and the algorithm makes no assumptions about types of phase transitions or 
trial wavefunctions. 

For our simulations we assume that the filaments are divided into an equal number of 
segments of equal length. This discretization leads to the Hamiltonian, 



H N (M) = Hf f (M) + H% lt (M) (A 1) 

where 

</ = a ^iM±lhM ( A2) 

j=i *;=i 

and 

M N N 

<-EEE^ lo s i ^ u) u) i , (A3) 

j=l k=l i>k 

and angular momentum 

M N 

J " = £Ew*tf)i 2 ' (a 4) 

3=1 k=l 

where 6 is the length of each segment, M is the number of segments, and Vfc(j) = 
x k{j) + Wk(j) is the position of the point at which two segments meet (called in PIMC 
a "bead") in the complex plane. 

The probability distribution for vortex filaments, 

G N (M) = , (A5) 



where 



Z N (M)= Gn(M), (A 6) 

allpaths 
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is the Gibbs canonical distribution. Our Monte Carlo simulations sample from this dis- 
tribution. 

The Monte Carlo simulation begins with a random distribution of filament end-points 
in a square of side 10, and there are two possible moves that the algorithm chooses at 
random. The first is to move a filament's end-points. A filament is chosen at random, 
and its end-points moved a uniform random distance. Then the energy of this new state, 
s' , is calculated and retained with probability 

A(s -> s') = min{l,exp (-[3[H™\s') ~ H%*(a)] - »[I N (s') - I N (s)})} , (A 7) 

where s is the previous state. (Self-induction, H s ^ 1 ^ , is unchanged for this type of move 
since it is internal to each filament.) The second move keeps end-points stationary and, 
following the bisection me thod of Ceperle y, grows a new internal configuration for a 



randomly chosen filament (jCeperlevi (|1995l )). Both the self-induction and the trapping 



potential are harmonic, so the Gibbs canonical distribution without interaction can be 
sampled directly as a Gaussian distribution. Therefore, in this move the configuration is 
generated by first sampling a free vortex filament, and then accepting the new state with 
probability 

A(s - s') = min {l,exp [-(3{H™\s') - H^\s))] } . (A 8) 

Our stopping criteria is graphical in that we ensure that the cumulative arithmetic 
mean of the energy settles to a constant. Typically, we run for 10 million moves or 50,000 
sweeps for 200 vortices. Afterwards, we collect data from about 200,000 moves to generate 
statistical information. 
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Appendix B. Spherical Model and the Saddle Point Method 



The spherical model was first proposed in a seminal paper of Berlin and Kac (jBerlin and Kac 



( 1952)), in which they were able to solve for the partition function of an Ising model given 
that the site spins satisfied a spherical constraint, meaning that the squares of the spins 
all added up to a fixed number. The method relies on what is known as the saddle point 
or steepest descent approximation method which is exact only for an infinite number of 
lattice sites. 

In general the steepest descent or saddle-point approximation applies to integrals of 
the form 

[ b e- N ^dx, (Bl) 

J a 

where f(x) is a twice-differentiable function, N is large, and a and b may be infinite. A 
special case, called Laplace's method, concerns real-valued f(x) with a finite minimum 
value. 

The intuition is that if xq is a point such that f(xo) < f(x)\fx ^ xq, i.e. it is a global 
minimum, then, if we multiply f{xo) by a number N, Nf(x) — N/(xq) will be larger 
than just f{x) — /(xo) for any x ^ xq. If N — > oo then the gap is infinite. For such large 
N, the only significant contribution to the integral comes from the value of the integrand 



at xq. Therefore, 



lim 



-Nf(x) 



1/N 



or 



( Berlin and Kad(|l952t ) 



lim log 

JV-^e N- 



dx = e-K X0 \ 



- Nf[x) dx = f(x ), 



(B2) 



(B3) 



Hartman and Weichmanl (|1995l )). A proof is easily obtained using 



a Taylor expansion of f(x) about xq to quadratic degree. 
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Appendix C. Evaluating the Free Energy Integral 



In this section we discuss our evaluation of the integral 



f[ir] = - log 



£>^exp(S) 



where 



5' 



/3'Llog(i? 2 )/4 



1 



daa'fS'l^^-] 2 + (it + 2[i)\iP(a)\ 2 - iR 2 r 
oo~ 



(CI) 



(C2) 



/?' = [3N , and a' = aN' 1 . 

The free-energy, Equation IC 1[ involves a simple harmonic oscillator with a constant 
external force, and we can re- write it, 



f[ir] = --irLR 2 - f3'L\og(R 2 )/4 - \nh\ir\. 



(C3) 



Here h is the partition function for a quantum harmonic oscillator in imaginary time, 



h[ir] — J D^exp ( J 



da-^m[\d^\ 2 +^\ 2 ]\ , 



(C4) 



which has the well-known solution for periodic paths in (2+l)-D where we have integrated 
the end-points over the whole plane as well, 



h[ir] = 



(e 



where m = a'f3' and cu 2 = [it + 2/i)/(a'/3') (jBrownl (|1992l) 



(C5) 



Zed boom 



Let us make a change of variables A = it + 2/i. Then the free-energy reads 
/[A] = (ji- \\)LR 2 - p'Llog(R 2 )/A - In 6 " 2 > 



(C6) 



where oj = y/\/{a'(3'). 
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